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S ! Abstract 

\ We present our results of lattice calculations of B parameters, which parameterize AB=2 

o ' 



transition amplitudes together with the leptonic decay constant. Calculations are made in the 
quenched approximation at (3=5.7, 5.9, 6.0 and 6.1, using NRQCD action for heavy quark and 
Q I the 0(a)-improved Wilson action for light quark. The operators are perturbatively renormal- 



ized including the correction of 0{as/{aM)"^) {m >0). We examine the scaling behavior of B 



^ ■ parameters, and discuss the systematic uncertainties based on the results with several different 

truncations of higher order terms in 1/M and Us expansions. We find Bs^imh) = 0.84(3)(5), 
BbJBb^ = 1.020(21)(1|^)(1^) and Bs,{mb) = 0.85(l)(5)(lo) in the quenched approximation. 
The errors represent statistical and systematic as well as the uncertainty in the determination of 
strange quark mass. 
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I. INTRODUCTION 



The determination of the Cabibbo-Kobayashi-Maskawa (CKM) matrix element \Vtd\ plays 
a crucial role in testing the unitarity relation of the CKM matrix, since the position of the 
vertex of the unitarity triangle would be essentially identified together with the angle 0i of 
the unitarity triangle. Now that the angle 0i has already been measured experimentally by 
the asymmetric B factories Q and its precision is expected to be improved substantially 
in the near future, the accuracy of \Vtd\ really determines the accuracy of the Standard 
model predictions. Then, other measurements of the CKM matrix elements, such as the 
determination of \Vub\ through a measurement of 6 — > ulu, may be used for a test of the 
CKM mechanism of the quark flavor mixing and the CP violation in the Standard Model. 

The CKM matrix element \Vtd\ may be determined using the mass difference AM^ in 
the neutral B meson mixing, as it emerges through a loop diagram mediated by top quark 
and W boson, which is proportional to iVf^^V^^p. The precision in the current world average 
(AMrf = 0.489 ± 0.008 ps^"*^ 0) is already as good as 1.6%. The constraint on \Vtd\ is, 
however, limited by the theoretical uncertainty in the calculation of the hadronic matrix 
element {B^\OlJB^) of the AB = 2 four-quark operator Ol^ = bj^{l — '^^)d b'y^{l — 
It is usually parameterized as ^f'^BsM'^ using the B meson leptonic decay constant fs and 
the B parameter Bb- In the vacuum saturation approximation, which is valid when both h 
and (anti-)(i quarks are nonrelativistic, Bb is normalized to unity. 

The best available theoretical method to calculate /s and -B_b is the numerical simulation 
of QCD on the lattice, whose current status is reviewed in ^ For the decay constant 
fBi several groups investigated the systematic errors in the lattice calculation, performing 
the simulations on several different lattices. It is found that the error associated with the 
large h quark mass is controlled reasonably well if one uses an effective theory for heavy 
quark such as the non-relativistic QCD (NRQCD) 0, |^ or the Fermilab formalism 0, and 



the results are insensitive to the lattice spacing |0 0, |T^, [T3|, |T4|. 

For the B parameter Bb-, on the other hand, most lattice calculations relying on the 
effective theory for heavy quark are limited to the static approximation, in which the l/m^ 
correction is neglected |jl5|, [1^, and the study of the systematic uncertainty depending 
on the lattice spacing has not been made. Recently some of us used the NRQCD action, for 
the first time, to calculate Bb 0], and the l/nib correction was studied at a fixed lattice 
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spacing |]T9|. They also calculated another B parameter Bs [|T9], ^ to parameterize the 
matrix element of the operator Os^ = b{l — 75)3 6(1 — 75)3, which appears in the heavy 
quark expansion of the width difference of B^ |21|, |22|, ^ . 

In this paper we extend the previous studies [jTB|, |T^, ^ to investigate the systematic 
errors in the calculation of the B meson B parameters. Using the same NRQCD action 
, |T^, ^ we calculate the B parameters at four lattice spacings to estimate the 



as m 



size of systematic errors depending on the lattice spacing. In order to minimize other 



sources of systematic errors, we use the 0(a)-improved Wilson quark action [24| for light 



quarks with the improvement coefficient csw calculated at the one-loop level ||2^, ^ ^ and 
nonperturbatively ^8 . 



Since NRQCD is an effective theory valid for heavy quark and the action is constructed by 
an expansion in inverse heavy quark mass, there is a potential source of systematic error due 
to the truncation of l/rrib expansion. Furthermore, in order to match the effective theory to 
the full theory one has to use perturbation theory, and errors from higher order corrections 
should also be taken into account. We introduce a method to estimate these systematic 
errors by treating the neglected higher order terms in different ways. It turned out that the 
error estimated in this way is quite consistent with a naive order counting assuming typical 
sizes for the expansion parameters. 

Due to the systematic errors discussed above, it is not straightforward to obtain an 
accuracy better than 10-15% for /bVBb which is relevant for the determination of \Vtd\- 
Alternatively, one could use the ratio AMg/ AM^, once the mass difference in the Bg — Bg 
mixing is measured. (The current experimental bound is AM^ > 13.1 ps^^ at 95% CL [p9(| .) 
It is proportional to i'^\VtslVtd\^ , where ^ is a ratio to describe SU(3) flavor breaking of the 



hadronic matrix elements given by 



Bs 



Since the bulk of the systematic errors in the calculations of and Bb cancels in this ratio, 
one may achieve much better accuracy, as stressed in |]30|. The largest remaining uncertainty 
comes from the chiral extrapolation of lattice data, which is also discussed in this paper. 

The B meson B parameters have also been calculated using the conventional relativistic 
actions for heavy quark ||31|, |3^. Since the lattice spacing in the present simulations is not 
small enough compared to the Compton wave length of the h quark, one has to extrapolate 
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the results obtained around the charm quark mass to the bottom quark mass, which is 
a significant source of systematic uncertainty. In fact, the extrapolation with the linear 
form in 1/M does not seem to agree with explicit calculations in the static limit p3| , |3^ , 
33. Therefore one may use the static result to constrain the heavy quark extrapolation 



in the infinite mass limit [BH]. We present a comparison of our result with these previous 
calculations. 

This paper is organized as follows. In the next section, we summarize some phenomeno- 
logical formulas for the mass and width difference in the — mixing. The lattice action 



and operators we employ in this work are defined in Section |T|, where the method to extract 
the continuum B parameters from lattice matrix elements is also described. Simulation de- 
tails and results are given in Section |I^ and 0, respectively. In Section ^ we present our 
results for the B parameters and their systematic uncertainties are discussed. Using these 
results we also predict the mass and width difference of B^ meson system. Calculation of 



the SU(3) breaking ratio ^ is briefly discussed in Section |VII| , and our conclusions are given 
in Section |VIII| . 



Preliminary reports of this work have already been presented in ||3^, |38|, p9|, HO 



II. PHENOMENOLOGICAL FORMULAS 

In this section we summarize the phenomenological formulas which involves the B meson 
B parameters. We also present some notations which will be used throughout the paper. 

A. Mass difference 



In the Standard Model the mass difference in the neutral B^ — B^ mesons (g denotes d 



or s) is given by 



mO,,{fi,)\B',). (2.1) 



So{xt) {xt = mf/m^) is the Inami-Lim function and ?72b is the short distance QCD 



correction [^, whose full expression is found in |^ 
The four-quark operator Oigifib) is deflned as 



OL^ilib) = &7m(1 - 75)g&7/.(l - 75)g, 



(2.2) 
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which depends on the renormahzation scale jj^b if it is defined in the continuum renor- 
malization schemes, such as the naive dimensional regularization (NDR) with the modi- 
fied minimal subtraction (MS) scheme. The scale dependence cancels with the prefactor 



7o/2/3o _|_ Q^/Jt) such that the physical mass difference is scale independent. In 
the NDR-MS scheme the anomalous dimensions are written as 

7o/9i 



and 



J n 



11 Uf, 

3 ^' 



71 



Pi 



102- -Uf, 
3 ^' 

4 



(2.3) 



(2.4) 



-7+ -Uf. 

9 ^ 



7o = 4, 7i 
The renormalization scale fib is usually taken at the b quark mass rub- 
The B parameter Bb„ is defined through 



{B'g\0,,ifibM) = -BsMfsMl, 



and the scale-independent ISb^ is given by 



Osifib) J 



(2.5) 



(2.6) 



The number of fiavors Uf is 5. To evaluate this expression we use the strong coupling 
constant as{fib) defined in the MS scheme with = 225 MeV, which corresponds to 
af\4.8 GeV) = 0.216. 

The bulk of the theoretical uncertainties cancels in the ratio 



AM, 

am: 



Vt 



ts 



td 



Bd 



(2.7) 



where ^ describes the SU(3) fiavor breaking of the matrix element {B^\0 Lq{fib)\Bq) as defined 
in ( |1 . 1|) . If one assumes the unitarity relation among the CKM matrix elements |VJs| — 
may be used to determine \Vtd\- 



B. Width difference 



Using the heavy quark expansion, the width difference in the neutral 5° — mixing is 
calculated as |^ ^ 

ATb, = -2^^{B,\ Imi I d^xTn,s{x)n,s{0)\B,), (2.8) 
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where TieS is the AB = 1 weak transition Hamiltonian. The main contribution comes from 
the transition bs cc followed by cc bs, and other contributions mediated by penguin 
operators are also considered. 

The operator product expansion (OPE) may be used to approximate the transition op- 
erator iJd'^xTT-Cc{i{x)T-Ccs{0), which gives an l/rrih expansion. At the leading order in l/rrib 
the AB = 2 four-quark operator Ol^ defined in (|2.2[) and another operator 



=6(1-75)5 5(1- 75)s, (2.9) 
appear. Then, the following formula 

^ XflA fl Mr, 

IK. 



/Ar\ ^ lQT:^B{B,^Xev)flMB 2 



\T J 9iz)r]QCD rri] 



'S^"^'""' ' ^"^-^^3 7^(mfe)2 



X [G{z)-BBXmb) + Gs{z)-^^^j^ + VT^z5y^\ (2.10) 



is obtained at the next-to-leading order [021, where mf, is the pole mass of b quark. The 



width difference is normalized by the total decay width of Bg meson F^^, which is written 
in terms of the semileptonic decay branching ratio B{Bs —>■ Xev) on the right hand side 
in order to remove an uncertainty in the value of \Vcb\- The phase space factor g{z) = 
1 — 8z + 8z^ — z'^ — Vlz^ Inz, where z = ml/ml, and the short distance QCD correction-*^ 

__ ^^2_31^^^_^^2 31 
37r 4 ^ ^ 2j ^ ^ 

are known factors, and the functions G{z) and Gs{z) describe the next-to-leading order 



2as{mb) 

Vqcd = 1 - 



QCD corrections ||2^ appearing in the calculation of the operator product expansion. The 



correction term 6i/m denotes the next-to-leading order contribution in the l/m^ expansion. 



which is estimated in using the factorization approximation. 



The B parameter Bs^ is defined through 



where 



{B'AOsMlB's) = -^fl^^j^^h (2-12) 

n{^^b) = ^^(^^)+^-(^^) (2.13) 
Mbs 

is the ratio of matrix elements of heavy-light axial vector current and pseudoscalar density 
and rh{fib) represents a MS quark mass. 



^ Following the treatment in the approximate form of |Q is used for fjqco- 
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In the following analysis, the scale jj^b is set to the pole mass of b quark, mi, = 4.8 GeV, for 
which G{z)=0.03 and Gs{z)=0.937. With input parameters z = 0.085, \Vcs\ = 1-^=0.976, 
Mb, = 5.37 GeV, 5(5, ^ Xeu) = 0.107, we obtain 



[ 0.0075B,(mb) + 0.2075s,(mb) - 0.077 ] . (2.14) 



\T Jbs V230 MeV^ 
For the central value of the decay constant, we choose a recent world average of unquenched 
lattice calculations fs, = 230(30) MeV 0. The uncertainties associated with these input 
parameters are discussed in Section |V1|. 



III. NRQCD ACTION AND OPERATORS 

In this section we describe the lattice NRQCD action and operators used in our cal- 
culations. The perturbative matching of the lattice operators to the continuum ones is 
summarized. 



A. NRQCD action 



To treat the heavy quark on a lattice with moderate lattice spacing a, the idea of the 



heavy quark effective theory (HQET) ^ |4^ is useful, as it allows us to describe the 
heavy quark of mass M without introducing large systematic errors scaling as a (positive) 
power of aM. In HQET the Lagrangian is organized as an expansion in inverse powers of 
M , and the terms beyond some fixed order are truncated. Since the physical expansion 
parameter is Aqcd/^ with Aqcd — 300-500 MeV, one may typically achieve the accuracy 
of order few percent for B mesons at the next-to-leading order, i.e. including terms of order 
1/M. 

Though NRQCD was originally introduced in the continuum ^ and on the lattice 
10, H to describe the quarkonium systems such as charmonium and bottomonium, for which 
the expansion parameter is a velocity of heavy quark rather than Kqqd/M, the formulation 



can also be used for the lattice study of heavy-light mesons as first demonstrated in |£0l . At 
the next-to- leading order in Aqcd/^, the Lagrangian in the continuum Euclidean space- 
time is written as 



rcont 
''-NRQCD 



Q 



Dn 



2M 



9- 



2M 



Dn 



2M 



9- 



a B 

2iW" 



X, 



(3.1) 
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for heavy quark field Q and heavy anti-quark field x- Both are represented by a two- 
component non-relativistic spinor. The derivatives Dq and D are temporal and spatial 
covariant derivatives respectively. The leading order term Dq represents a heavy quark as a 
static color source. The leading correction terms of order Aqcd/^ are the non-relativistic 
kinetic term D^/2M and the spin- (chro mo) magnetic interaction term a ■ B/2M, where B 
denotes the chromomagnetic field strength. In the usual HQET approach, only the leading 
terms are kept in the effective Lagrangian and corrections of order Aqcd/jW^ are treated as 
operator insertions. Alternatively, in our lattice calculation we include the correction terms 
in the Lagrangian (|3.1| ) and evaluate the matrix elements with the heavy quark propagator 
including the effect of order Kqqj^/M. 

On the lattice we use a discretized version of the Lagrangian whose explicit form 

is written as 



•^NRQCD = Y.Q^(^)^^^,y - KQ{x,y))Q{y) + Y.x\x){6^,y - K^{x.y))x{y). 



(3.2) 



x,y 



x,y 



The kernel to describe the time evolution of (anti-)heavy quark is given by 



KQ{x,y) 



1 - 



aHn 



1 



2n 
aHo\" 



1 



a6H^ 



a6H' 



1 - 



aHn 



2n 



{x,y), (3.3) 
{x,y), (3.4) 



2n J \ 2 J ^ \ 2yV 2n 
where n denotes a stabilization parameter introduced in order to remove the instability 
arising from unphysical momentum modes in the evolution equation [0, H. The operator 
^4^^ is defined as 6'^\x,y) = Sxi±i,yiSx,y, and Hq and 6H are lattice Hamiltonians defined 
by 

A(2) 



^^0 = 

5H ^ 



2aMo' 
9 

'^^2aMa 



-a B, 



(2) 



(3.5) 
(3.6) 



where A^^^ = J2f=i ^1 is the Laplacian defined on the lattice with A 
and A^""* being forward and backward covariant derivatives in the i-th direction. In ( p.6|) 
the chromomagnetic field operator B is the usual clover-leaf type lattice field strength . In 
these definitions, the lattice operators A*^^^ and B are dimensionless, i.e. appropriate powers 
of a are understood. The space-time indices x and y are implicit in these expressions. The 
bare heavy quark mass Mq is distinguished from the renormalized one M. 
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At the tree level, the lattice action ( p.2|) describes the continuum NRQCD in 
the limit of vanishing lattice spacing a. (We may identify Mq with M and take the tree 
level value cb = 1-) The leading discretization error for the spatial derivative is of order 
{akQQ-oY Kqqd / M . Since the temporal derivative is discretized asymmetrically, the leading 
error appears at order aD^, whose typical size is estimated as oAq^j-j/M^ using the equation 
of motion. The gauge potential part is automatically improved, as it is exponentiated into 
the temporal link variable f/4. 

In the presence of radiative corrections, the heavy quark mass Mq and the chromomag- 
netic coupling cb have to be tuned in such a way that the continuum values are reproduced 
at each value of the strong coupling constant a^. Furthermore, the radiative corrections 
generate many other terms which do not exist in the continuum Lagrangian ( |3.1|) , because 
NRQCD is not a renormalizable field theory. In general these terms appear with some factor 
of form Q;J/(aMo)'" with positive integers k and m {k>l and m > 2).^ Therefore, NRQCD 
should be considered as an effective theory valid for small 1 / (aM) up to higher order terms 
in l/(aM). 

Perturbation theory can be used to calculate the renormalization of the parameters. For 
example, the one-loop calculations of energy shift and mass renormalization were carried out 
for lattice NRQCD by Davies and Thacker |^ and by Morningstar and by ourselves 
|TB|, |TP[ for the above particular form of the NRQCD action.^ To improve the perturbative 



expansion we utilize the tadpole improvement namely, all the gauge links in the action 
( p. 21) are divided by its mean field value uq determined from the plaquette expectation value 
Mo = ((Trf/p)/3)^''^, where some counter terms are introduced in perturbative calculations. 
The one- loop tuning of the coupling constant Cb for the spin-(chromo)magnetic interaction 
term ( p.6|) has not yet been performed, so we take the tree level value cb = 1 after making 
the tadpole improvement. 

The relativistic four-component Dirac spinor field h is related to the two-component non- 
relativistic field Q and x appearing in the NRQCD action ( |3.2| ) via the Foldy-Wouthuysen- 



^ There are also the lowest dimension operators Q^Q and x^X; but they only give the energy shift and do 
not contribute to the dynamics of heavy quark. 



3 



We note that the evolution kernels (|3.3| ) and (3^) are slightly different from the definition used, for 



example, in |53], where the (1 — aHo/2n)^ terms appear inside of the (1 — aSH/2) terms. 
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Tani (FWT) transformation 

' 7 . A(±) 



1 - ^ h, (3.7) 

2aMo / ^ ' 



where 

' Q 



(3.8) 

The symbol A*^^) denotes a symmetric covariant differentiation operator Aj-^'' = A^-'*''' A-~\ 



B. Bilinear operators 

The heavy-light axial-vector current Ap, = b^y^'y^q and pseudoscalar density V = h^^^q 
appear in the definition of the B parameters through the vacuum saturation approximation. 
We use the calligraphic symbols and V to denote the currents defined in the continuum 
full theory. Since the pseudoscalar density diverges in the continuum, it is renormalized 
with the MS scheme at a scale /i. On the other hand, the axial-vector current does not need 
renormalization, because it is partially conserved in the continuum full theory. 

The corresponding lattice operators are 

jP = iTq, (3.9) 

^'' = 2jiH^-^'^')r«. (3.10) 

where Jr is ^4 for T = 7475 or P for F = 75. The light quark field q is described by the 
0(a)-improved Wilson quark action [^]. We apply the tadpole improvement for the 
light quark field using the critical hopping parameter Kc to define the mean link variable 



Mo = 1/8kc, so that we normalize the light quark field with a factor Jl — j^. The heavy 



quark field h is defined in 

The one-loop matching between the continuum and lattice operators is written as 



Jt = 

with one-loop coefficients pp"*. The coefficient p^^ is written as 



4°^ + ^P?4n^p?)4^ (3.12) 



= 2ln{a'M') + CA, (3.13) 
pP = I Mp^/M') + ^ Ma'M^) + Cp, (3.14) 
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for r = 7475 and 75, respectively. In the static limit the numerical constants are (a = —16.55 



51, m, m and Cp = -11.21 EM. For the NRQCD action (O) with a finite heavy quark 



mass Mq the numerical values for (^a and are available in Table III of |T9|. 

In the static limit, while the second term of ( |3.12| ) vanishes, the third term remains 
finite and describes the ©(asaAqco) improvement, and its coefficient pp'' is p""^ /2aMQ = 
Pp /2aMQ = 13.01 [Q. Away from the static limit, these terms give contributions of the 
©(asaAqco) and 0{asAQCT)/M), and the one-loop coefficients are calculated only for the 
axial vector current F = 7475 for our choice of the NRQCD action T^.^ 



C. AB=2 operators 

We assume that the continuum four-quark operators Ol{p) and Os{p) are renormalized 
in the MS scheme with totally anti-commuting 75. In the renormalization of Osip), the 
subtraction of evanescent operators is made with the definition given by Eqs. (13)-(15) of 



22 



For matching of these four-quark operators, the following lattice operators are involved 
at the lowest dimension: 

Ol = 67M(l-75)g&7M(l-75)g, (3.15) 

Or = h^{l + -f,)q h^{l + -f,)q, (3.16) 

Os = 6(1 -75)g&(l- 75)9, (3.17) 

= 2 67^(1 - 75)g 67^(1 + 75)g + 4 6(1 - 75)5 6(1 + 75)?, (3.18) 

Om = 2 67^(1 - 75)g 67m(1 + 75)g - 4 6(1 - 75)9 6(1 + 75)?, (3.19) 

Op = 2 67^(1 - 75)g 67^(1 + 75)g + 12 6(1 - 75)? 6(1 + 75)5, (3.20) 

Ot = 5 67^(1 - 75)g 67^(1 + 75)g - 34 6(1 - 75)5 6(1 + 75)5. (3.21) 

As in the bilinear operator case, one has to introduce higher dimensional operators to remove 
0(a) errors at the one- loop level. For the four-quark operators, however, the necessary one- 
loop calculations to remove the ©(a^aAqcD) error are made only in the static limit |2^, |4 



^ The same quantity was previously calculated in ^ , but for a slightly different NRQCD action. 

^ Note that a different notation is used in [T^ . Similar calculation was previously made by Morningstar 

and Shigemitsu 
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The higher dimensional operators appearing in that hmit are 



Old = &7/^(l - Ts)? &7m(1 " 75)(a7 ■ A(=^))g, 
Ond = 2 67^(1 -75)g 67^(1 + 75)(a7-AW)g 

+4 6(l-75)g6(l + 75)(«7- AW)g, 
O5D = 6(l-75)g&(l-75)(a7- AW)g, 
Opi5 = 2 67^(1 - 75)g 67/.(l + 75)(a7 ■ A(^))g 

+12 6(l-75)g6(l + 75)(a7- AW)g. 



(3.22) 

(3.23) 
(3.24) 

(3.25) 



The one-loop matching is written as follows: 



Ol + ^Pl,lOl + ^Cl,sOs 



+-;—Cl,rOr + —C,l,nOn + --;—Ql,mOm 
Att Att Atx 

+ -r-CL,LDOLD + -T—Ql.NdOnD, 

An An ' 

Os + T^Ps.sOs + -t^Ps^lOl 
An An 

+-;—(s,rOr + -T^Cs,pOp + -t^Cs^tOt 
An An An 



a. 



a. 



a. 



+ ^C,s,sdOsd + --;—Cs,ldOld + --;—Cs,pdOpd, 
An An An 



(3.26) 



(3.27) 



where the coefficients Pl,l, Ps,s, and ps,L contain the physical scale /i and as follows. 



Pl,l 
Ps,s 

PS,L 



-2\nl^]+A\n{a'M') + CL,L, 



3 J 



(3.28) 
(3.29) 
(3.30) 



The numerical results for the one-loop coefficients (l,x {X = L, S, R, N, and M) and (s,x 
{X = S, L, R, P, and T) are given in Table VI and Vlll of |jT9[. In the static limit, Cl,m 



and Cs,T vanish, and others agree with the previous calculations |T^, 0, ^ ^ 01 

The last hues in ( p.26| ) and ( p.27| ) are added to remove the error of O^asaAqcri), but 
their coefficients are known only in the static limit. Their values are Cl,ld = —17.20, 
Cl,nd = -9.20, Cs,sd = -6.88, Cs,ld = 2.58 and Cs,pd = 1.15 |2D|, g. 



A numerical error in was later corrected in |5^, 61 
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D. Truncation of expansions 



As in the matching of the NRQCD action discussed in Section [III A| , we have to truncate 
the 1/M and the perturbative expansions in the matching of the bihnear and four-quark 
operators. The 1/M expansion is truncated at 0(1/M), which is consistent with our choice 
of the NRQCD action (pTTP, and the perturbative corrections of order and higher are 
neglected. 

In addition, there are mixed corrections of 0{asAQCD/M). In the matching of the bi- 
linear operators the matching coefficients for the mixed corrections are available and such 



corrections were actually included in [O by combining with higher dimensional operators as 



shown in ( p.l2|) . For the four-quark operators, however, the mixing with higher dimensional 
operators at the one-loop level has not been calculated yet''. Thus, in this paper, the mixed 
corrections are not considered in both of the bilinear and four-quark operators. This means 
that, for the bilinear operators, only the first term of (|3.12|) is taken, thus the matching 



becomes multiplicative in this approximation. 

At this level of accuracy, it is arbitrary to apply the FWT transformation to a heavy 
quark field, if the heavy quark field forms an operator appearing in the one-loop corrections. 
Namely, in ( p. 261) and (|3.27|) , we may replace all the four-quark operators Ox except for 
that in the first term by O'x, where 

O'l = h,{l - l5)q - l,)q. (3.31) 

and so on, and the heavy quark field h is not rotated by the FWT transformation ( p.Tj) . 
Therefore O'x differ from Ox at 0(Aqcd/^)- In the naive order counting analysis both 
choices are equivalent up to unknown corrections of 0(q;sAqcd/A^^)- 
In the calculation of B parameters through the ratios 

{B\OM\B) 

Bs{f^>^) - |(5|P(^,)|0)(0|P(;.,)|5)' ^'-''^ 

the perturbative and 1/M expansions may be truncated in several different ways. A natural 
choice to match the ratios (|3.32|) and ( |3.33|) is to write the numerator and denominator as 



^ Except for the static limit, where the mixing terms describe the correction of order agQ rather than ag/M. 
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they stand: 



1 + I^PL,L 



jDlat I Os /■ jDlat 



B,X 



1 I Os „ nlat _i_ as „ TDlat 1 «s /■ TDlat 

^ + i^PS.S + 4^PS,L-Ds^L + 2^X=R,P,T 4^^S,X^s,X 



1 + 

J- + 4^Pp 



(3.34) 
(3.35) 



A roman numeral, (I) in this CclSG, clS db superscript oi Bb oi Bs distinguishes the method to 
truncate the expansion. -B^*j^ and B^^^-^ are B parameters defined with the lattice operators 
as 



B 



lat 
B,X 



{B\Ox\B) 



B 



lat 
S,X 



l{B\Af^\0){0\AP\By 
{B\Ox\B) 



(3.36) 
(3.37) 



|(5|P(o)|0)(0|P(o)|5)' 

which are directly measured in the numerical simulation from a ratio of correlation functions 
as we describe in the next section. 

Alternatively, one may linearize the perturbative expansion as 



bT\p. 



l + ^(p.,.-2pr) 
l + ^{ps,s-2pP) 



B^L + 



E 

X=S,R,N,M 



Air 



/■ Tjlat 
^L,X^B,Xi 



(3.38) 



Bts + ^Ps,LB's:L+ E £C5,x5g^. (3.39) 

X=R,P,T 



Formally they are different from the method I by order a^, which is not known. We ex- 
pect, however, that perturbative expansion behaves better for the method II, because the 
contributions from factorized diagrams to the four-quark operators are the same as those of 
the corresponding bilinear operators, so that the radiative corrections partly cancel in the 
combination p^,^ - 2pf or ps,s - 2pg\ 

For each method I or II, we also consider the variation of replacing Ox and J^^ by O'x 
and J'P , respectively, as discussed above, and define the methods as I' and IF, where J'^ 
are defined similarly to ( |3.31| ). The difference of the method F (IF) from I (II) is of order 
a.AqcD/M. 

Since the level of accuracy of these four methods is equivalent in the naive order counting 
argument, they can be used to estimate possible systematic errors due to the truncation of 
expansions. 
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IV. LATTICE SIMULATIONS 



A. Simulation sets 

We have performed numerical simulations at four j3 values. For three of them {(5=6.1, 
5.9, and 5.7), which we call the simulation set A, the 0(a)-improvement coefficient csw 
in the light quark action is determined using the one-loop expression csw = + 
0.199av'(l/'^)]- The one-loop coefficient is calculated in and we apply the tadpole 



improvement with the plaquette expectation value to define the mean link variable. For 
the last simulation (/3=6.0), which we denote as the simulation set B, the non-perturbative 
value is used for csw HH- Therefore, as far as the light quark sector is concerned, the 
discretization error is minimized in the set B, for which the leading error is of O(a^), while 
the effect of O(a^a) is remaining in the set A. For the quantities involving heavy quarks, 
however, both sets of simulations give the same order of accuracy, since the heavy quark 
action and operators are not improved at the same level. 

Simulation parameters are summarized in Table |. The parameters for the simulation 



set A are almost the same as in our previous work for the leptonic decay constant |]13 
except that the number of statistical ensembles is increased in this work to obtain stable 
signals for three-point functions. The set B is our new simulation set, which is intended for 
comparison with our recent unquenched simulations [^, and its (3 value, /3 = 6.0, is 
chosen so that the inverse lattice spacing becomes about 2 GeV. In this paper, we present 
only the quenched results leaving the unquenched calculations for future publications. 

For both simulation sets, A and B, we take the standard plaquette gauge action, and 
the configuration generation and gauge fixing are made as in [Q. Four values of the light 
quark hopping parameter are chosen for each (3 as given in Table |. They correspond to the 
light quark mass rUg covering the range 771^/2 < rriq < 2ms, where denotes the physical 
strange quark mass. The hopping parameter corresponding to the strange quark mass is 
determined using the K or (p meson masses as input, and will be denoted as Hgi and k,s2 
respectively. 

The heavy quark mass in our simulation ranges from 2mb/3 to Arrib- The smallest heavy 
quark mass in the lattice unit aMo is limited around unity due to the possibly large system- 



atic error in the matching calculation as discussed in Section III. The limit in the heaviest 
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side is set by the exponentially growing statistical error . 

The lattice spacing a is determined through the string tension (for the set A) or the rho 
meson mass (for the set B). For the simulation set A it is confirmed that both determinations 
are in good agreement (3.5% variation depending on (3) [0. Therefore, in effect the lattice 
spacing is set using the rho meson mass for both data sets. 

We use the simulation set B to obtain our central value and the other to investigate 
the systematic errors depending on the lattice spacing. The primary reason is that the 
discretization error is minimized by the non-perturbative 0(a) improvement for the set B. 
The set B is also advantageous since we have larger statistics and hence the numerical results 
are more stable. 



B. Correlation functions 



The method to calculate two- and three-point functions mostly follows that of [|T^. We 
put a local source at the origin of the lattice and solve for the light quark propagator. The 
heavy quark and anti-quark propagators are obtained from the same local source by solving 
the evolution kernels (|3.3| ) and (|3.4| ), respectively. 

Three point functions are constructed as 



C?Hii,i2) =EE(0|^ (ti,fi)Ox(0,0)Af (t2,X2) |0) 



l5t 



(4.1) 



where Ox is one of the four-quark operators defined in ( p.l5|) -( p.21D . We take ti > and 
^2 < so that a B meson propagates in the positive direction in time and a B meson 
propagates in the opposite direction. We also measure two point functions 

Cfit) = J2{0\r[Afit,x)Af\0,0)] |0), 
cPit) = E(0|T[Af(t,x)p(°)(0,0)] |0), 

X 

for positive and negative values of t. 
A smeared current Af, defined as 



(4.2) 
(4.3) 



^4{t, ^) = E Hy)Kt, X + y)7475g(^, x), 



(4.4) 



is used to enhance the overlap with the ground state B meson. We measure the smearing 
function 0(x) for each set of heavy and light quark masses with a limited number of gauge 
configurations before starting the main simulation. 
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We extract the lattice B parameters -B^x ( P-36|) and B^^x ( P-37|) from the following ratios 



for large enough \ti 



RB,xiti,t2) = 



Cx\ti,t2) 



,(3) 



B 



lat 
B,Xi 



-^Cp {ti)Cp (t2) 



1,2). 



(4.5) 
(4.6) 



C. Meson masses 

In order to calculate the heavy-light meson masses precisely, we also calculate two-point 
functions with the smeared source and local sink, 



Cf'-'it) = ^(0|T[Ar(t,f)Af(0,0)] |0), 

X 

cP'-'it) = J2m\Af^\t,x)P'{0M |0), 



(4.7) 
(4.8) 



for which the statistical signal is much better than those with the local source and smeared 
sink. The heavy-light meson mass is, then, obtained by adding the binding energy -Ebin mea- 
sured from the two-point functions and bare quark mass oMq. Including one- loop corrections 
we use 

aMp = Z^aMo + Eun - Sm, (4.9) 
where perturbative corrections Zm and 6m are given as 



Zm = 1 + asB, 



6m = a^A, 



(4.10) 
(4.11) 



and A and B for each bare quark mass are given in Table. I of [13 



V. SIMULATION RESULTS 



A. Ratio of correlation functions 



We first extract the B parameters defined on the lattice Bp^x ^^s^xi which are 
obtained from the asymptotic behavior of the ratios -R_b,x(^i, ^2) and Rs,x(ti,t2) as (f4.5|) - 
([4.6|). In Figures (top and middle panels in each Figure) we plot these ratios as a 
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function of ti for some fixed values of ^2- For illustration we show the operators giving 
leading contributions, i.e. -B^*^ and -B^^g, for the heavy quark mass closest to the physical b 
quark mass and the lightest quark mass (largest k, value). We obtain similar plots for other 
mass parameters, but the statistical signal becomes much noisier for larger heavy quark 
mass. 

The range of t {ti and ^2) to be included in the fit of the ratios has to be chosen such 
that the effect of excited states is negligible. We identify the plateau seen in the plots of 
RB,x(ti, ^2) and Rs,x(ti, ^2) as the region where the ground state contribution dominates. To 
be more conservative, we also check that the plateau is reached for the effective mass plot 
of two-point functions C^\t) and C^^^^{t), which are calculated for the same smearing 
function as used in the calculation of three-point functions. The plots are shown in the 
bottom panel of Figures 

In the fit of the ratios we take a range of t as wide as possible in order to avoid possible 



contamination from the statistical fluctuation |6J]. The fit is done for a fixed value of 
t2 = hmin and changing ti in the range [timin,timax]- The value of [timin,timax] is hsted in 
Table |. In order to quantify the possible effect from excited state contamination, we also 
carried out a fit with larger values of timm (= ^2)- Since the statistical error grows rapidly 
as timin is taken larger, the maximum change for timin is chosen to keep the statistical error 
smaller than 8-10%. The effect for Bb is found to be 1% or less except for (3 = 6.1 where it 
is at most 3%. For Bs it brings a 1.5-3% effect except for (3=5.7 where it is negligible. The 
variations among different choices of the fit range are taken into account in the final results. 



B. Chiral extrapolation 

The results of 5^*^ and 5if^ are insensitive to the light quark mass. An example is 
shown in Figure ^ where the data at (3 = 6.0 are plotted as a function of amg = 1 (^^ — -^^ . 
The B parameters for all relevant operators are plotted: X = L, R, S, N and M for -B^X) 
and X = S, L, R, P, T for -B^^f- Each B parameter is normalized by its vacuum saturation 
approximation. We averaged the matrix elements with X = L and R, as they should be 
equal in the infinite number of statistics by parity symmetry. In later sections, the averaged 
matrix elements are denoted by X = LR. 

We extrapolate these B parameters to the chiral limit of light quark assuming a linear 
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function in amq. In most cases the chiral extrapolation changes the value of B parameters 
from the lightest measured data by about 1% or less. Therefore, the chiral extrapolation is 
extremely stable and the associated systematic error is negligible. To confirm this observa- 
tion we also tried a quadratic extrapolation for some parameter sets, for which we find that 
the results are consistent with the linear extrapolation within the statistical error. 

In chiral perturbation theory for heavy-light mesons, the logarithmic dependence such as 
rriq Inrriq is predicted for Bb ||65|, In the quenched approximation the chiral limit is even 
divergent as Inm^. The more divergent term Inm^ has a factor 1 — 3g^ as its coefficient, and 
the B*B7i coupling g is evaluated in the range 0.2-0.7 using D* Dn decay |6^, D* decay 
width |6^, quark models and quenched lattice calculations |j70|, |7I|, |7^. It means that 
this divergent logarithm is relatively unimportant because of its small coefficient 1 — Sg"^ = 
0.2(7). It is, however, difficult to resolve such logarithmic dependences from the data taken 
in the range of our light quark masses. In this work, therefore, we do not further consider 
them, leaving the study of the chiral logarithm including the effect of unquenching for future 
pubhcations. 

Results of the linear extrapolation are summarized in Table [n|-pi^, where we list the 
values of -B^x and B''gx at each (3 and qMq. The value of k corresponding to the physical u 
or d quark mass, which we denote Kud, is very close to the critical value Kc- The value of Kg 
corresponding to the strange quark mass depends on the input quantity. We list the results 
at Ksi, for which the K meson mass is used as input, and at k,s2, for which meson mass is 
used. 



C. 1/Mp dependence 

The 1/Mp dependence of the lattice B parameters -B^x (-^sx) is plotted in Figure ^ 
(Figure |^). The light quark is extrapolated to the chiral limit. Although the data at 
different (3 values are overlaid, they do not necessarily agree because the operators are not 
renormalized. Comparison with an adequate definition (the B parameter in the continuum 



renormalization scheme) will be discussed in detail in Section [VF . 

We find that the mass dependence is small for -B^lj?) while it is significant for others. This 
behavior can be mostly understood using the vacuum saturation approximation (VSA) ||18| 



In VSA the matrix element in the numerator of the B parameter is generally decomposed 



19 



into |(0|A^|P)p and |(0|P|P)p. For 5^^^,^, however, it is written by |(0|y4^|P)p only and no 
I (0|P|P) p term appears by definition, so that B^^^^ = 1 is independent of 1/Mp. For others, 
the term | (0|P|P) ^ gives a strong mass dependence proportional to {Mp/MY = (1+A/M)^, 
where A represents the binding energy produced by the light degrees of freedom. Comparison 
of the lattice data with VSA is made in Ref. |TB| . 



D. Renormalized B parameters 

The B parameters for the continuum operators are obtained from the lattice B param- 
eters using ( |3.34| ) for Bb and ( p.35| ) for Bs- We consider the truncation method I in this 



subsection. The results of other truncations are discussed in the next subsection. 
In order to see the effect of 1/M corrections we consider the quantity 

^BA^^^)^{^^}]"^^ Bs{^^,) (5.1) 
\as[MB) ) 

as a function of 1/Mp. The factor {Q.s{Mp)las[MB)f'^^'^ is introduced to cancel the loga- 
rithmic dependence on M coming from the continuum one-loop integral, so that the heavy 
quark expansion in 1/Mp is explicit. Up to two- loop corrections $B^(/ib) is equivalent to 
Bsil^'b) obtained with a replacement of M in p^ L ( p.28|) and in ( p.l3| ) by the physical 
b quark mass nib, which can be confirmed by expanding the factor (as{Mp) /as{MB)y^'^° in 
as{Mp) explicitly. With the replacement the static limit (1/Mp —>■ 0) of NRQCD simply 
becomes the conventional static approximation. Therefore, in the calculation of ^BBif^b) we 
explicitly set the physical b quark mass rrib in the matching coefficients ( |3.28[ ) and ( p.l3| ). It 



should be also noted that at the physical B meson mass, namely Mp=Mb, our definition of 
^Bsil^h) exactly agrees with the definition (p. I]). 

Figure |] (top panel) shows ^BBif^b) at /3 = 6.0. The light quark mass is extrapolated 
to the chiral limit, and the renormalization scale pb is set to rUb. In the one-loop matching 
( p. 34 ) we use the renormalized coupling ay(g*) defined through the heavy quark potential 



|. The scale q* represents the momentum region where the relevant one-loop integral 
dominates. Since it is not known, we use three typical values it /a, 2/ a and 1/a and consider 
their variation as an indication of systematic uncertainty from higher order perturbative cor- 
rections. We find that the variation among different coupling constants becomes substantial 
as one goes to the static limit, while it is relatively unimportant in the physical mass region 
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1/Mp ~ 0.2 GeV ^. This is because the one-loop coefficients in the matching ( p.34|) grows 



toward the static hmit. For the 1/Mp dependence of ^BBif^b) we observe a shght positive 
slope and curvature, but the large systematic uncertainty implies that the mass dependence 
is insignificant. 

We obtain a similar plot for ^Bsif^t) in Figure ^ (bottom panel), which is an analog of 
^BgifJ'b) but for Bs{fib)- The definition of ^Bsif^b) with the renormalization group improve- 
ment as in { f).l\j is more complicated, since the logarithmic dependence appears in more than 
one coefficients, so that we have to consider a mixing of operators. In this work, however, 
we avoid this problem by replacing the heavy quark mass in ps^s and ps,L by the physical 
value as we did for (^Bgil^h)- 

For Bsifib) the one- loop coefficients are relatively small and their dependence on the heavy 
quark mass is mild. Hence, we obtain a smaller variation due to different scale settings in 
the coupling constant. 

E. Effect of truncation of expansions 



As we discussed in Section [I11D| there are several method to truncate the perturbative 
and 1/M expansions. We consider the following four methods. In the methods I and I' 
the perturbative matching is truncated in the numerator and denominator separately as in 



( p.34| ) and ( p.35|) , while in the methods II and IF the denominator is linearized as (|3.38|) 



and ( p.39|) . In the primed methods the heavy quark field without the FWT rotation 
is used for one-loop correction terms. 

In Figure | we plot (^Bgil^b) (top panel) and (^Bsil^b) (bottom panel) for four different 
truncation methods. As can be seen from the figure, the methods I and I' (or II and IF) 
agree in the static limit, since their difference is only in the FWT rotation. On the other 
hand, the difference between the methods I and II (or F and IF) is smaller for lighter heavy 
quarks because the one-loop correction in the denominator becomes small. 

We consider the variation among different truncations as an indication of systematic 
uncertainties from higher orders of perturbative and 1/M expansions. The error estimation 
is given in the next section. 
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F. Results at different lattice spacings 



Since NRQCD is formulated by an expansion in 1/M and not a renormalizable field 
theory, it does not allow a continuum limit; instead it has to be considered as an effective 
theory valid in the region where l/(aMo) is small enough. The dependence of systematic 
errors on the lattice spacing a is not just a simple power series in a, but contains its inverse 
powers. Therefore, the question is how one can find a region of a where the discretization 
error is small while the errors scaling as 1/a" (n is a positive integer) is under control. 
Although the order counting argument as discussed in []TB|, ^ provides a rough estimate of 
errors, it is essential to confirm it using actual simulation data. 



In Figure |T0| we plot ^Bsif^b) (top panel) and ^Bsif^t) (bottom panel) obtained with the 
truncation method 1 for four different lattice spacings. The largest (smallest) inverse lattice 
spacing is 2.3 GeV at P = 6.1 (1.1 GeV at /? = 5.7). We find that around the physical B 
meson mass {1/Mb ~ 0.2 GeV^^) the results agree within order 10% for ^BBif^b) or even 
better for $Bs(/^fe)- The agreement becomes marginal toward the static limit especially for 
the nonperturbatively improved lattice (/5 = 6.0), but it is not statistically significant. 

Results of physical ^^(mb) (top panel) and Bs{mb) (bottom panel) are plotted in Fig. ^ 
as a function of the lattice spacing, where the variation due to the different choice of fit range 
is added to the error bar at each jS. The a dependence is a mixture of the discretization 
error scaling as a"^ and the truncation error containing a form like 1/a". In addition, 
the truncation of perturbative expansion gives a functional dependence like 1/lna. It is, 
therefore, difficult to determine the shape of the a dependence, but the data imply that none 
of these errors is diverging in the region we measured. 



VI. PHYSICS RESULTS 



A. Analysis of systematic errors 

As discussed in the previous section, we have performed the calculation of the B parame- 
ters with four different truncations of 1/M and ag expansions. Furthermore, the calculations 
are made at four different lattice spacings. All of these calculations have different amount 
of various systematic errors, and thus they allow us to estimate the uncertainty in our final 
results. In this subsection we first list possible sources of systematic errors and estimate 
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their size using a naive order counting. Then, their results are compared with the actual 
lattice data. 

One of the possible systematic errors arises from the discretization of derivatives, which 
scales as a power of the lattice spacing a. Because our actions and operators are 0(a)- 
improved at tree level, the leading error is of order a? and of order a^a. Since we are 
using an effective theory for heavy quark, the truncation of the 1/M expansion leads to a 
systematic error. For our choice of actions and operators the leading contribution is of order 
Again, since the matching of the 1/M terms is done at the tree level only, we also 
expect an error of order a^/Mfc. The perturbative matching of operators are truncated at 
the one-loop level, so that there is an uncertainty of order a]. In Table ^ we list their typical 
size at each (3 value using a naive order counting. Where the scale is needed we assume the 
typical spatial momentum inside a heavy-light meson to be p ~ Aqcd ~ 300 MeV. For the 
strong coupling constant we use a typical value av'(2/a) as listed in Table |. 

The contribution of order /M^ is not investigated in this paper, as we just neglect 
the terms in the action and operators. In Table ^we assign 2% as the corresponding 

uncertainty rather than a naive order counting 0.4%, taking the estimate from explicit lattice 



study in |18 



Since we have removed errors of as/ (aM)"^ (m >0) by perturbative matching, the leading 
contribution which prevents us from the continuum extrapolation with the NRQCD action 
has the form a1/{aM). Although its size in the naive order counting is smaller than the 
pure two-loop correction a^, we include it in our error analysis (thus in Table as it gives 
the leading contribution growing toward the continuum limit. 



As mentioned in Section [III D| , the results from the different truncation methods (I, II, I' 
and ir) are expected to differ from each other by 0{a1) or 0{asp/Mb). We compare their 
results in Figure |1^ for Bb and in Figure |13| for Bs- The results of the four truncation 
methods when as{2/a) is used in the one- loop matching are plotted. In these Figures we 
also show the size of the systematic errors of order 0(0;^) and 0{asp/Mb) estimated with 
the naive order counting (first two lines of Table ^ added in quadrature). Although the 
statistical error of the data points makes the comparison somewhat ambiguous, we conclude 
that the naive order counting reasonably explains the scatter of the lattice data. The same 
conclusion is reached when we use different values of q* in the perturbative matching. 

An alternative way to estimate the 0{al) error is to see the variation of results with 
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q*=Ti /a^ 2/ a and 1/a. From Figure H, where the data at /3 = 6.0 are plotted, we find that 
the variation with the choice of scale in the coupling constant is consistent with our order 
counting (~ 3%) for the data points around the B meson mass. In the static limit, on the 
other hand, the one-loop coefficient is uncomfortably large, and the variation among the 
results with different q* is much larger than our naive estimate. 

Other sources of systematic errors to be tested are O(a^p^), 0{apas) and 0(a^/(aM{,)). 
Although it is difficult to disentangle various a dependent systematic errors solely from 
the data, our results are stable against the change of lattice spacing suggesting that the 
systematic error is well estimated by the naive order counting. 

Finally we also investigated the systematic error from the contamination of excited states 
by changing the fit range for the ratio RB,x(ti,t2) and Rs,x(ti,t2)- We find that the effect 
is at most 3%. In particular, at (3=6.0 it is found to be 1% for Bb and 3% for Bs- These 
variations are taken into account in the final results. 



B. B parameter results 

As we discussed above, the systematic errors estimated with the naive order counting 
actually describe the differences among different calculations. We, therefore, use the order 
counting argument to quote our estimate of systematic uncertainties. 

We take the central values from the data at /5 = 6.0 (set B) with q*=2/a, and obtain the 
following results in the quenched approximation. 

BBim,) = 0.84(3)(5), (6.1) 

BBAmb) = 0.86(2)(5)(0), (6.2) 

Bsirub) = 0.82(2)(5), (6.3) 

BsAmb) = 0.85(1)(5)(1^), (6.4) 

where the first and second errors represent statistical and systematic ones respectively. In 
the systematic error the contamination from excited states is added in quadrature with 
other sources estimated by the order counting. The third error is from the uncertainty of 
rris arising from the different input physical quantities, i.e. itlk or m^. 

The corresponding renormalization scale independent B parameter is obtained from ( p.6|) 



Bb = 1.29(5)(8), (6.5) 
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which may be compared with the previous calculations using the relativistic actions for 
heavy quark, 1.38(11)(1!]) ^ and 1.40(5)(i?) [gg. These two resuhs are slightly higher 



than our result, although is still consistent within the large error. A possible reason for 
the high values in the relativistic approach is in the extrapolation in the heavy quark mass 
from the charm quark mass region to the bottom. In fact, the combined analysis of the 
HQET and relativistic heavy quark, in which the interpolation in 1/M can be made, gives 
1.34(6) (t^) m, z.e. closer to our result. 



C. Applications 

In this subsection we present a few examples of physics applications of our results. It 
should be noted, however, that our calculation is still in the quenched approximation and 
there is no rigorous estimate for the associated uncertainty. For the following applications 
we assume that the quenching effect is negligible for the B parameters as suggested by our 



preliminary calculations [pHl E^, M 



For the B meson decay constant, on the other hand, the large effect of quenching has 
been found @, ^, ^j. Furthermore, a large uncertainty due to the presence of chiral logarithm 



is suggested for fs |73|, while the effect is not too large for fs^. In the following analysis 
we therefore consider the quantities for which only Jb^ is needed, and use the recent world 
average of unquenched lattice calculations of fs, = 230(30) MeV @, ^ when needed. 

Assuming the three generation unitarity \Vts\ — we obtain the mass difference in 

5° — 5° mixing using ( p.l|) as 

AM, = 19.4(5.5) ps-\ (6.6) 

where the statistical and systematic errors in theoretical and experimental quantities are 
added in quadrature, but the final error is dominated by the uncertainty of /b^- The value 
is consistent with the current lower bound AM, > 13.1 ps~^ at a 95% CL |2^. Tevatron 
Run II is expected to measure the mass difference very precisely in a few years. 

The width difference in the Bs meson system could also be measured at Tevatron Run 



II, if it is large enough. Using (|2.14|) and ( |6.4D , we obtain 



/AT 



) = 0.106 ±0.020 ±0.028 ±0.037 ±0.024, (6.7) 
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where the first through third errors are from m?, = 4.8(3) GeV, fs^ = 230(30) MeV, and TZ 
(or Thb{mb) = 4.25(25) GeV and rhs{mb) = 0.10(3) GeV), respectively. The last error comes 
from the uncertainty in the estimation of the 1/m correction, for which we assign 30%. The 
error from the B parameters is much smaller than the others listed above. 

The uncertainty in the calculation of (Ar/r)^^ is still very large (~ 50% if added in 
quadrature). In order to improve it one has to calculate the 1/m corrections reliably, as it 
largely cancels the leading contribution from Bs as seen in (|2.14|) . Currently, only an upper 
bound is obtained for this quantity experimentally. Our prediction is consistent with the 
bound (^)ij, < 0.31 at a 95 % CL ||. 



VII. SU(3) BREAKING RATIO ^ 

Since we expect that the bulk of systematic uncertainty in the calculation of fs and Bb 
cancels in their SU(3) breaking ratio fss/ fs and BbJBb, they could be useful to reduce 
the errors in the determination of \Vtd\ through the relation ( |2.7| ). 

In the lattice calculation, the deviation of their ratio from unity is the quantity to be 
calculated and the errors scale as BbJBb — 1 rather than BbJBb itself. In the present 
case, the naive estimate of SU(3) breaking is 0(m|-/A^) ~ 25%, where A^ is a scale of the 
chiral symmetry breaking ~ 1 GeV, and the order counting of uncertainties for the ratio are 
starting from this order. 

As done in Section [VI A| , we compare our order counting with data in Figure |l^. As 



we expected, the variation with different truncations is much smaller than the expected 
systematic errors (dotted line). The dependence on (3 is sizable, but not significant compared 
to the relatively large statistical error. 
Our result is 

:^ = 1.020(21)(i}^)(^^), (7.1) 

where the central value is taken from (3 = 6.0 (set B), and the errors are statistical, systematic 
and the uncertainty of in the order given. 

For the calculation of ^ ( |1.1|) we have to combine BbJBb with fBs/fs, which is called 
the indirect method. On the other hand, it is also possible to directly obtain ^ from a ratio 
of {OlJ and {OlJ (direct method). It has been discussed that one may obtain rather large 
value of ^ if one uses the direct method BO, R^. Therefore, in the following we check if we 
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could obtain consistent results from both methods using our data. 



Figure |T5| shows the chiral extrapolation of (O^), as required in the direct method. The 
data are obtained at /? = 6.0 (set B) for a heavy quark mass closest to the b quark mass. 
The dashed line is obtained by a linear fit to the data (open circles), while the solid curve 
represents a fit with a linear plus quadratic term in anig. Although the data look consistent 
with the linear fit, the chiral limit with the quadratic fit is higher by about one standard 
deviation. 

An open diamond at the chiral limit, on the other hand, is obtained through the indirect 
method, i.e. the decay constant and B parameter are separately extrapolated to the chiral 
limit with a linear fit. Although we have not presented a calculation of the decay constant 
in this paper, they are done on the same set of gauge configurations ai (3 = 6.0 and the 



lattice axial current is renormalized as described in Section [IIIB . The result is completely 



consistent with the quadratic fit in the direct method. It implies that in the direct method 
the chiral extrapolation is more difficult and needs enough statistics to control, since (Ol) 
is effectively the decay constant squared so that the finite aniq correction is amplified^. 

VIII. CONCLUSIONS 

In this work, we calculate the B meson B parameters on the lattice in the quenched 



approximation. The calculation is an extension of our previous works ||T8| , pO| , in which Bb 
and Bs were calculated for the first time with the lattice NRQCD action. 

In the present work we include a detailed study of systematic uncertainties. Using the 
lattice simulations at four different (3 values with the 0(a)-improved actions, we find that the 
B parameters are essentially insensitive to the discretization error. We also investigate the 
systematic errors associated with the truncation of heavy quark and perturbative expansions, 
which are necessary in the effective theory approach such as NRQCD. By comparing four 
different truncations of these expansions, we are able to confirm that the naive order counting 
argument of the systematic errors could actually give a reasonable estimate. 

In our final results for the B parameters the systematic error is ~ 6%, which is already 
smaller than that in the equivalent calculations of fs (10-20%), owing to the fact that it is 



The similar discussion may be found in [p2[ 
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defined as a ratio to the vacuum saturation approximation. Further reduction of systematic 
errors, if it is necessary, requires higher order calculation of perturbation theory and the 
O(a^) improvement. Approaching to the continuum limit will not help to reduce the errors 
in the NRQCD approach. 

For a precise extraction of the important CKM element \Vtd\ through the SU(3) breaking 
ratio of the B — B mixing one needs It is preferable to take the chiral limit separately 
for Jb and Bb, as they have milder light quark mass dependence. The SU(3) breaking ratio 
of is obtained with accuracy of order a few percent, since the B parameter is extremely 
insensitive to the light quark mass and the large cancellation of systematic errors is expected. 

The largest remaining uncertainty in our calculation is in the quenching approximation, 
though it is not explicitly discussed in the paper. We are currently performing an unquenched 
simulation with the same lattice action at similar lattice spacing, which will allow us to 
directly study the quenching effect. 
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set 




B 


/3 


6.1 


5.9 


5.7 


6.0 


size 


X 64 16^ X 48 12^ x 32 


20^ X 48 


#conf 


518 


419 


420 


655 


csw 


1.525 


1.580 


1.674 


1.769 


1/a (GeV) 


2.29 


1.64 


1.08 


1.82 


K 


13586 


13630 


13690 


0.13260 




0.13642 


0.13711 


0.13760 


0.13331 




0.13684 


0.13769 


0.13840 


0.13384 




0.13716 


0.13816 


0.13920 


0.13432 




0.13635 


0.13702 


0.13800 


0.13355 




0.13609 


0.13657 


0.13707 


0.13318 




0.13767 


0.13901 


0.14157 


0.13531 


(aMo,n) 


(7.0,2) 


(10.0,2) 


(12.0,2) 


(10.0,2) 




(3.5,2) 


(5.0,2) 


(6.5,2) 


(5.0,2) 




(2.1,2) 


(3.0,2) 


(4.5,2) 


(3.0,2) 




(1.5,3) 


(2.1,3) 


(3.8,2) 


(2.1,3) 




(0.9,4) 


(1.3,3) 


(3.0,2) 


(1.3,3) 




[8,261 


[6,171 


[4,131 


[ 7,18] for aMo = 10.0 
[ 9,18] for aMo = 5.0 
[10,18] for aMo = 3.0, 2.1, 1.3 


Mo 


0.8816 


0.8734 


0.86087 


0.87603 


ay(7r/a) 


0.149 


0.164 


0.188 


0.159 


0^/(2/0) 


0.172 


0.193 


0.229 


0.186 




0.229 


0.270 


0.355 


0.256 



TABLE I: Simulation parameters. For the simulation 
is determined at the tadpole-improved one-loop level, 
value [28] is used. 



set the 0(a)-improvement coefficient csw 
For the set i?, the nonperturbatively tuned 
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aMo 


12.0 6.5 4.5 3.8 3.0 


Tjlat 

J-' DTD 

Kud 
Ksl 
Ks2 


0.916(21) 0.906(17) 0.894(15) 0.889(14) 0.882(14) 
0.931(13) 0.923(11) 0.915(09) 0.911(09) 0.905(09) 
0.934(12) 0.927(09) 0.920(08) 0.916(08) 0.911(08) 


TDlat 

Kud 
Ksl 
Ks2 


-0.656(11) -0.708(10) -0.765(10) -0.802(11) -0.870(12) 
-0.659(07) -0.713(06) -0.770(07) -0.806(07) -0.872(08) 
-0.659(06) -0.714(06) -0.772(06) -0.808(06) -0.872(08) 


jylat 

JD D AT 

Kud 
Ksl 
Ks2 


1.220(36) 1.459(31) 1.707(32) 1.864(34) 2.144(40) 
1.212(23) 1.444(20) 1.683(21) 1.831(23) 2.095(27) 
1.210(20) 1.440(17) 1.676(19) 1.822(20) 2.081(24) 


TDlat 

JD o Byf 

l^ud 
Ksl 
Ks2 


-6.44(11) -6.87(10) -7.30(09) -7.58(10) -8.08(10) 
-6.44(07) -6.87(06) -7.31(06) -7.58(07) -8.06(07) 
-6.44(07) -6.87(06) -7.31(06) -7.58(06) -8.06(07) 



TABLE II: Numerical values for lattice B parameters /3=5.7. 
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aMo 


12.0 6.5 4.5 3.8 3.0 


JDo a 
Kud 
Ksl 
Ks2 


0.931(15) 0.909(11) 0.894(10) 0.887(10) 0.879(09) 
0.941(10) 0.927(07) 0.916(07) 0.910(06) 0.904(06) 
0.944(09) 0.931(07) 0.922(06) 0.917(06) 0.910(05) 


r>lat 

Kud 
Ksl 
Ks2 


-1.301(31) -1.163(23) -1.044(19) -0.982(18) -0.891(17) 
-1.330(19) -1.200(14) -1.088(12) -1.028(12) -0.938(11) 
-1.338(17) -1.209(13) -1.099(11) -1.040(10) -0.951(10) 


B¥p 

Kud 
Ksl 
Ks2 


-12.59(18) -12.56(15) -12.52(14) -12.51(13) -12.51(12) 
-12.66(12) -12.68(10) -12.69(09) -12.69(09) -12.70(08) 
-12.68(10) -12.71(09) -12.73(08) -12.74(08) -12.76(08) 


Kud 
Ksl 
Ks2 


55.43(87) 54.09(67) 52.91(58) 52.31(55) 51.45(50) 
55.79(57) 54.72(45) 53.77(39) 53.27(37) 52.50(35) 
55.89(51) 54.90(40) 54.01(36) 53.53(34) 52.78(32) 



TABLE III: Numerical values for lattice B parameters Bg-^ at /3=5.7. 
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aMo 10.0 

r>lat 
^B,LR 

Kud 0.977(61) 

Ksi 0.956(34) 

Ks2 0.952(30) 



Dlat 
I^s2 



TDlat 
^B,M 

Ks2 



-0.669(32) 
-0.650(17) 
-0.646(15) 



-6.715(297) 
-6.521(176) 
-6.476(157) 



5.0 



0.936(35) 
0.931(19) 
0.930(16) 



-0.718(19) 
-0.715(10) 
-0.714(09) 



1.595(74) 
1.560(41) 
1.552(35) 



-7.156(182) 
-7.087(104) 
-7.071(091) 



3.0 



0.904(25) 
0.911(13) 
0.913(11) 



-0.800(16) 
-0.805(09) 
-0.806(08) 



1.979(60) 
1.935(35) 
1.924(31) 



-7.903(162) 
-7.862(093) 
-7.852(082) 



2.1 



0.884(22) 
0.897(12) 
0.899(10) 



-0.904(17) 
-0.911(10) 
-0.912(09) 



2.453(67) 
2.381(41) 
2.364(37) 



-8.828(177) 
-8.742(102) 
-8.723(090) 



1.3 



0.848(22) 
0.871(12) 
0.877(11) 



-1.148(25) 
-1.150(17) 
-1.150(15) 



3.541(102) 
3.395(066) 
3.361(059) 



-10.904(241) 
-10.684(145) 
-10.634(128) 



TDlat 
^B,N 

Hud 1.346(122) 

Ksi 1.288(070) 

Ks2 1.275(060) 



TABLE IV: Numerical values for lattice B parameters 5^*^ at /3=5.9. 
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aMo 


10.0 5.0 3.0 2.1 1.3 


JDo Q 

Kud 
Ksl 
Ks2 


0.942(44) 0.907(23) 0.882(15) 0.868(13) 0.848(11) 
0.923(24) 0.912(12) 0.899(08) 0.888(07) 0.873(06) 
0.918(21) 0.913(11) 0.902(07) 0.893(06) 0.878(06) 


r>lat 

Kud 
Ksl 
Ks2 


-1.376(85) -1.184(45) -0.999(30) -0.850(25) -0.627(22) 
-1.357(48) -1.188(24) -1.018(16) -0.875(15) -0.662(14) 
-1.353(42) -1.188(21) -1.022(14) -0.880(13) -0.670(12) 


b^t^ 

Kud 
Ksl 
Ks2 


-13.24(55) -13.08(31) -13.10(24) -13.20(22) -13.33(22) 
-12.90(31) -13.02(17) -13.10(13) -13.17(12) -13.27(12) 
-12.82(27) -13.01(15) -13.10(12) -13.17(11) -13.25(10) 


b,l 

Kud 
Ksl 
Ks2 


57.69(233) 55.84(126) 54.53(96) 53.71(89) 52.33(88) 
56.36(138) 55.71(072) 54.75(55) 53.86(50) 52.36(49) 
56.05(122) 55.68(063) 54.80(48) 53.90(43) 52.36(42) 



TABLE V: Numerical values for lattice B parameters Bg^ at /3=5.9. 
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aMo 7.0 

TDlat 
^B,LR 

Kud 0.833(66) 

Ksi 0.892(37) 

Ks2 0.904(33) 



Ks2 



TDlat 
Ks2 



-0.653(30) 
-0.658(18) 
-0.659(16) 



-7.125(338) 
-6.857(220) 
-6.804(204) 



3.5 



0.846(32) 
0.892(18) 
0.901(16) 



-0.701(17) 
-0.725(10) 
-0.730(09) 



1.674(61) 
1.689(37) 
1.692(34) 



-7.690(180) 
-7.575(112) 
-7.552(102) 



2.1 



0.845(23) 
0.882(12) 
0.889(11) 



-0.808(14) 
-0.836(09) 
-0.842(08) 



2.179(55) 
2.181(35) 
2.181(33) 



-8.545(143) 
-8.452(089) 
-8.433(082) 



1.5 



0.834(22) 
0.871(12) 
0.878(10) 



-0.929(16) 
-0.958(11) 
-0.964(10) 



2.715(66) 
2.700(44) 
2.697(41) 



-9.490(155) 
-9.389(099) 
-9.369(092) 



0.9 



0.819(24) 
0.857(13) 
0.865(12) 



-1.255(28) 
-1.282(20) 
-1.287(19) 



4.130(115) 
4.047(078) 
4.031(073) 



-12.046(230) 
-11.877(156) 
-11.843(146) 



TDlat 
^B,N 

Hud 1.227(112) 

Ksi 1.294(064) 

Ks2 1.307(057) 



TABLE VI: Numerical values for lattice B parameters 5^*^ at /3=6.1. 
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Tjlat 



«s2 



7.0 



3.5 



2.1 



1.5 



0.9 



0.892(40) 0.843(20) 0.834(13) 0.831(12) 0.827(11) 
0.907(24) 0.879(11) 0.869(7) 0.863(7) 0.855(7) 
0.910(21) 0.886(10) 0.876(7) 0.870(6) 0.861(6) 



r>lat 



Ks2 



-1.138(92) -1.016(41) -0.870(26) -0.741(22) -0.536(19) 
-1.229(52) -1.080(23) -0.916(15) -0.784(13) -0.571(12) 
-1.247(46) -1.094(21) -0.925(14) -0.792(12) -0.578(12) 



TDlat 
^S,P 



-13.09(54) -13.27(26) -13.31(19) -13.32(18) -13.35(17) 
-13.01(35) -13.28(17) -13.31(12) -13.31(11) -13.31(11) 
-12.99(33) -13.28(15) -13.31(11) -13.31(10) -13.31(10) 



Ksl 
I^s2 



58.65(260) 56.91(124) 55.23(85) 53.90(76) 51.70(68) 
57.32(171) 56.66(078) 55.11(53) 53.78(47) 51.65(42) 
57.06(158) 56.61(072) 55.08(48) 53.76(43) 51.63(38) 



TABLE VII: Numerical values for lattice B parameters Bg-^ at (3=6.1. 
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aMo 


10.0 5.0 3.0 2.1 1.3 


Tjlat 

J-' DTD 

Kud 
Ksl 
Ks2 


0.820(53) 0.832(51) 0.868(39) 0.857(30) 0.858(25) 
0.864(37) 0.869(29) 0.888(22) 0.877(17) 0.872(15) 
0.874(34) 0.877(25) 0.892(18) 0.881(14) 0.875(13) 


TDlat 

Kud 
Ksl 
Ks2 


-0.574(28) -0.649(27) -0.738(23) -0.820(20) -1.011(22) 
-0.601(21) -0.673(17) -0.757(14) -0.840(12) -1.029(14) 
-0.607(19) -0.678(15) -0.761(12) -0.845(11) -1.033(13) 


jylat 

JD D AT 

Kud 
Ksl 
Ks2 


1.128(107) 1.376(97) 1.787(84) 2.207(75) 3.095(81) 
1.148(073) 1.393(57) 1.785(48) 2.184(45) 3.023(53) 
1.152(066) 1.397(49) 1.784(41) 2.178(39) 3.007(48) 


TDlat 

JD o Byf 

l^ud 
Ksl 
Ks2 


-6.67(36) -7.33(28) -7.76(24) -8.44(20) -10.04(20) 
-6.68(24) -7.15(17) -7.69(14) -8.35(12) - 9.88(13) 
-6.68(23) -7.11(15) -7.67(12) -8.33(11) - 9.84(12) 



TABLE VIII: Numerical values for lattice B parameters -B^x /3=6.0. 
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aMo 


10.0 5.0 3.0 2.1 1.3 


JDo a 
Kud 
Ksl 
Ks2 


0.829(40) 0.847(35) 0.854(26) 0.840(19) 0.830(14) 
0.868(30) 0.881(22) 0.881(15) 0.869(11) 0.857(08) 
0.877(28) 0.888(20) 0.887(13) 0.875(10) 0.863(07) 


r>lat 

Kud 
Ksl 
Ks2 


-1.184(77) -1.085(68) -1.005(48) -0.879(33) -0.705(23) 
-1.249(53) -1.138(39) -1.034(27) -0.907(19) -0.726(14) 
-1.264(50) -1.149(34) -1.040(23) -0.913(16) -0.730(13) 


B¥p 

Kud 
Ksl 
Ks2 


-12.87(58) -13.17(47) -13.13(37) -13.17(28) -13.32(22) 
-12.95(41) -13.02(27) -13.11(22) -13.14(17) -13.26(14) 
-12.97(38) -12.99(24) -13.11(19) -13.14(15) -13.24(12) 


Kud 
Ksl 
Ks2 


57.82(276) 58.01(211) 55.59(161) 54.30(118) 52.93(88) 
57.99(194) 56.99(126) 55.44(094) 54.21(072) 52.78(57) 
58.02(180) 56.76(111) 55.41(082) 54.19(063) 52.74(51) 



TABLE IX: Numerical values for lattice B parameters -B^^ at /3=6.0. 
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/3 


5.7 


5.9 6.1 6.0 




5% 


4% 3% 3% 


0{asp/M^) 


1% 


1% 1% 1% 




8% 


3% 2% 3% 


0{apas) 


6% 


4% 2% 3% 


0{al/{aMb)) 


1% 


1% 1% 1% 






2 % 


Total (added in quadrature) 


11 % 


7% 5 % 6 % 



TABLE X: Estimate of systematic uncertainties by a naive dimensional counting at each /? value. 
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FIG. 1: RB,L{ti,t2) (top) and Rs,s{tiM) (middle) at /3 = 5.7, oMq = 3.8 and k = 0.13920. 
Horizontal line represents a fit with a range ti = [4, 13] for a fixed t2 = 4. The bottom plot shows 
an effective mass for two point functions C^^\t) (open symbols) and C^^^^{t) (filled symbols). 



Circles and squares represent data points for positive and negative t respectively. 
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FIG. 2: Same as Fig. ||, but for /? = 5.9, oMq = 2.1 and k = 0.13816. Horizontal line represents 
a fit with a range ti = [6, 20] for a fixed t2 = 6. 
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FIG. 4: Same as Fig. ||, but for /? = 6.0, uMq = 2.1 and k = 0.13432. Horizontal line represents 
a fit with a range ti = [10, 18] for a fixed t2 = 10. 
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FIG. 5: Chiral extrapolation of S^*^ (top panel) and .B^x (bottom panel) at /3 = 6.0 and uMq 



lat 



2.1. Data are normalized by their vacuum saturation approximation (VSA). It is S^]^'^^ = 1, 
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B,S 



-5/8, B 
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B,N 



1, 5 
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B,M 



-6 for 5^*x, and Bgg 



(VSA) 



1 r(VSA) 



-8/5, = 



64/5, Sg"^^^ = 288/5 for In VSA the correction of order 1/M is neglected. 



-ylat 



47 




1/M„ (GeV" ) 



0.1 0.2 

1/M„ (GeV) 



FIG. 6: 1/Mp dependence of the lattice B parameters Bb,x {X = LR, S, N and M). A quadratic 
fit is plotted for the data at (3 = 6.0. 
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FIG. 7: 1/Mp dependence of the lattice B parameters Bs,x {X = S, LR, P and T). A quadratic 
fit is plotted for the data at (3 = 6.0. 
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FIG. 8: ^Bgil^'b) (top panel) and ^Bsil^b) (bottom panel) for fij, = mi,. Data show the result at 
P = 6.0 in the chiral limit of the light quark. The truncation method I is chosen as a demonstra- 
tion. Different symbols correspond to different scales of the coupling constant in the perturbative 
matching, and data points are slightly shifted in the x-direction for clarity. The data at the static 
limit (1/Mp = 0) is obtained by an extrapolation in 1/Mp with a quadratic function. 
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FIG. 9: ^Bgil^b) (top panel) and ^Bsif^b) (bottom panel) for fij, = m^. Data show the result 
at /? = 6.0 in the chiral limit for the light quark and renormalized with the coupling Q;y(2/a). 
Different symbols correspond to different truncations of perturbative and heavy quark expansions. 
The data at the static limit (1/Mp = 0) is obtained by an extrapolation in 1/Mp with a quadratic 
function. 
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FIG. 10: ^Bgil^b) (top panel) and ^Bsif^b) (bottom panel) for fij, = nib. Results at different 
lattice spacings are compared. Data show the results in the chiral limit for the light quark and 
renormalized with the coupling ay(2/a). The truncation method I is chosen as a demonstration. 
The data at the static limit (1/Mp = 0) is obtained by an extrapolation in 1/Mp with a quadratic 
function. 
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FIG. 11: Dependence of BBirrih) (top panel) and Bs{mi,) (bottom panel) on the lattice spacing 
a. Data show the results in the chiral limit for the light quark and renormalized with the coupling 
ay(2/a). The truncation method I is chosen as a demonstration. The variation due to the different 
choice of fit range is added to the error bar at each (3. 
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FIG. 12: Results for Bb at each (3 with four different truncation methods (see text). Comparison 
is made with the estimate using the naive order counting (a band given by the dotted hnes). 
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FIG. 13: Same as Figure |l^, but for Bs- 
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FIG. 14: Same as Figure but for BbJBb 
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FIG. 15: Comparison of direct and indirect methods. The dashed hne is a hnear fit and the solid 
curve is obtained with a fit with linear and quadratic terms. An open diamond at the chiral hmit 
is obtained through the indirect method. 
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